######
#This portion of the code simply merges ideal point and alliance data onto our events data#
######

rm(list=ls())
setwd("~/Dropbox/Foreign Gifts/Replication materials")
library(devtools)
#The version of record of this paper uses the countrycode package as of the git commit dated April 2, 2015 and found at https://github.com/vincentarelbundock/countrycode/commit/893d2319f542e411ab70f117ae5d25709747012e  Earlier versions of the countrycode package contain errors that are consequential for merging the data.  It is conceivable that later versions will also impact merging behavior, although we believe that this version is accurate for the purposes of the operations we conduct.  Nonetheless, replicators should install that version from github to reproduce the results here.
if (F) {install_github('vincentarelbundock/countrycode',ref="893d2319f542e411ab70f117ae5d25709747012e")}
library(countrycode)

## load('internationalEvents.RData')
internationalEvents <- read.csv("internationalEvents.csv")

#We use Version 7.1 of the Bailey, Strezhnev, and Voeten dyadic ideal point data (date May 11, 2015).  Replicators may independently obtain this file at https://dataverse.harvard.edu/dataset.xhtml?persistentId=hdl:1902.1/12379&version=7.1  This version of the data must be used  to reproduce the results here.  The file loaded here is an unmodified version of the original.

## load('Idealpointsdyadicdistance_1.RData')
x <- read.csv("Idealpointsdyadicdistance_1.csv")

nonUSData<-internationalEvents[internationalEvents[,"actor1country"]!="UNITED STATES"&internationalEvents[,"actor2country"]!="UNITED STATES",]
nonUSData[,"actor1country"]<-countrycode(nonUSData[,"actor1country"],"country.name","country.name")
nonUSData[,"actor2country"]<-countrycode(nonUSData[,"actor2country"],"country.name","country.name")

idealPoints<-x[x[,"countryname1"]=="United States of America",c("year","cowcode2","idealpointdistance")]
idealPoints<-data.frame(idealPoints,"my"=floor((idealPoints[,"year"]-1)/4))
aggDist<-aggregate(idealpointdistance~my+cowcode2,data=idealPoints,FUN=mean)
colnames(aggDist)[3]<-"aggDist"
idealPoints<-merge(idealPoints,aggDist)
idealPoints<-idealPoints[,c("year","cowcode2","aggDist")]
colnames(idealPoints)<-c("year","actor1country","idealpointdistance.1")
idealPoints[,"actor1country"]<-countrycode(idealPoints[,"actor1country"],"cown","country.name")
nonUSData<-merge(nonUSData,idealPoints,all.x=T)
colnames(idealPoints)<-c("year","actor2country","idealpointdistance.2")
nonUSData<-merge(nonUSData,idealPoints,all.x=T)

#This is an unmodified copy of version 4.1 of the COW alliances dataset.
alliances<-read.csv("COWAlliances.csv")
#Note that the US appears first in all alliances of which it is a member in this data
usalliances<-alliances[(alliances[,"ccode1"]==2&alliances[,"year"]>1945),]
usalliances[,"state_name2"]<-countrycode(usalliances[,"ccode2"],"cown","country.name")
usalliances<-usalliances[,c("state_name2","year","defense")]
usalliances<-aggregate(defense~state_name2+year,data=usalliances,FUN=sum)
colnames(usalliances)<-c("actor1country","year","ally.1")

#Note to replicators: This merge results in a "1" for ally.1 and ally.2 whenever a country IS allied with the United States, but often produces an "NA" when a country is NOT allied with the United States.  An alternative merge should be used to run an analysis on countries NOT allied with the US (no such results are presented in the paper or supplementary information)
nonUSData<-merge(nonUSData,usalliances,all.x=T)
colnames(usalliances)<-c("actor2country","year","ally.2")
nonUSData<-merge(nonUSData,usalliances,all.x=T)

nonUSData$ally.1[is.na(nonUSData$ally.1)] <- 0
nonUSData$ally.2[is.na(nonUSData$ally.2)] <- 0

nonUSData <- subset(nonUSData, !is.na(idealpointdistance.1) & !is.na(idealpointdistance.2))

library(countrycode)
library(parallel)
library(lubridate)
library(xts)
library(nlme)
library(lmtest)
library(tseries)
library(sandwich)
## library(stargazer)


cutoff<-median(c(nonUSData[,"idealpointdistance.2"],nonUSData[,"idealpointdistance.1"]))
inUSBloc<-nonUSData[nonUSData[,"idealpointdistance.1"]<=cutoff&nonUSData[,"idealpointdistance.2"]<=cutoff,]
#
usBlocTS<-xts(inUSBloc[,"GOLDSTEINSCALE"],
              order.by=as.Date(paste(inUSBloc[,"year"], substr(inUSBloc$date, 6, 7), "01", sep = "-"))
                  )
monthlyMean<-apply.monthly(usBlocTS,mean)
monthlyCount<-apply.monthly(usBlocTS,nrow)
isElection<-year(monthlyCount)%%4==0
inWindow<-month(monthlyCount)%in%c(6,7,8,9,10,11)
monthlyCount<-data.frame(monthlyCount,isElection)
proximateElection<-isElection&inWindow
proximateMidterm<-year(monthlyMean)%%4==2&inWindow
theMonth<-month(monthlyMean)
theYear<-year(monthlyMean)
augmented<-cbind("meanScale"=monthlyMean,"proximate"=proximateElection,"premidterm"=proximateMidterm,"Month"=theMonth,"Year"=theYear)
reelect<-augmented[,"proximate"]&(!augmented[,"Year"]%in%c(1952,1960,1968,1988))
norelect<-augmented[,"proximate"]&(augmented[,"Year"]%in%c(1952,1960,1968,1988))

augmented<-cbind(augmented,NA)
#Fixed effects are assembled at the President-term level - i.e., switch to a new FE
#When President changes (incl. via death) or after a subsequent innauguration
colnames(augmented)[ncol(augmented)]<-"Term"
augmented["1933-03-01/1937-01-31","Term"]<-1
augmented["1937-02-01/1941-01-31","Term"]<-2
augmented["1941-02-01/1945-03-31","Term"]<-3
augmented["1945-04-01/1949-01-31","Term"]<-4
augmented["1949-02-01/1953-01-31","Term"]<-5
augmented["1953-02-01/1957-01-31","Term"]<-6
augmented["1957-02-01/1961-01-31","Term"]<-7
augmented["1961-02-01/1963-11-31","Term"]<-8
augmented["1963-12-01/1965-01-31","Term"]<-9
augmented["1965-02-01/1969-01-31","Term"]<-10
augmented["1969-02-01/1973-01-31","Term"]<-11
augmented["1973-02-01/1974-08-01","Term"]<-12
augmented["1974-08-01/1977-01-31","Term"]<-13
augmented["1977-02-01/1981-01-31","Term"]<-14
augmented["1981-02-01/1985-01-31","Term"]<-15
augmented["1985-02-01/1989-01-31","Term"]<-16
augmented["1989-02-01/1993-01-31","Term"]<-17


# Table 11
model4<-lm(meanScale~reelect+norelect+premidterm+as.character(Month),data=augmented)
coeftest(model4,vcov = vcovHAC(model4))
#
model6<-lm(meanScale~reelect+norelect+premidterm+as.character(Month)+as.character(floor((Year-1)/4)),data=augmented)
coeftest(model6,vcov = vcovHAC(model6))
#
model5<-lm(meanScale~reelect+norelect+premidterm+as.character(Month)+as.character(ceiling(Year/2)),data=augmented)
coeftest(model5,vcov = vcovHAC(model5))

## stargazer(
##     coeftest(model4, vcov=vcovHAC(model4)),
##     coeftest(model6, vcov=vcovHAC(model6)),
##     coeftest(model5, vcov=vcovHAC(model5)),
##     report=('vcsp'),
##     type="latex",
##     omit=c("as.character","Constant"),
##     digits=2
##     )
## nobs(model4); nobs(model6); nobs(model5)

presidentConflict <- aggregate(meanScale ~ Term, augmented, FUN = mean)
names(presidentConflict)[2] <- "AvgConflict"
augmented <- merge(data.frame(augmented), presidentConflict, by.x="Term",by.y="Term")

model.interaction<-lm(scale(meanScale/AvgConflict)~proximate*scale(Year)+as.character(Month),data=augmented)
coeftest(model.interaction,vcov = vcovHAC(model.interaction))




#BLS Recessions
augmented<-cbind(augmented,F)
colnames(augmented)[ncol(augmented)]<-"Recession"
augmented["1937-05-01/1938-06-01","Recession"]<-TRUE
augmented["1945-02-01/1938-10-01","Recession"]<-TRUE
augmented["1948-11-01/1949-10-01","Recession"]<-TRUE
augmented["1953-07-01/1954-05-01","Recession"]<-TRUE
augmented["1957-08-01/1958-05-01","Recession"]<-TRUE
augmented["1960-04-01/1961-02-01","Recession"]<-TRUE
augmented["1969-12-01/1970-11-01","Recession"]<-TRUE
augmented["1973-12-01/1975-03-01","Recession"]<-TRUE
augmented["1980-01-01/1980-07-01","Recession"]<-TRUE
augmented["1981-07-01/1982-11-01","Recession"]<-TRUE
augmented["1990-07-01/1991-03-01","Recession"]<-TRUE

FEMA <- read.csv("FEMA_AllYears.csv")
FEMA$Date <- as.Date(FEMA$Date, "%m/%d/%Y")

FEMA <- xts(
    FEMA,
    order.by=as.Date(paste(year(FEMA$Date), month(FEMA$Date), "01", sep = "-"))
    )
FEMA <- apply.monthly(FEMA, nrow)
names(FEMA) <- "disasters"

augmented <- merge(
    augmented,
    FEMA
    )

augmented$disasters[is.na(augmented$disasters)&index(augmented)>="1953-01-01"] <- 0

## Table 17
model1<-lm(meanScale~proximate+Recession++as.character(Month)+as.character(floor((Year-1)/4)),data=augmented)
coeftest(model1,vcov = vcovHAC(model1))
model2<-lm(meanScale~proximate+scale(disasters)+as.character(Month)+as.character(floor((Year-1)/4)),data=augmented)
coeftest(model2,vcov = vcovHAC(model2))

## stargazer(
##     coeftest(model1, vcov=vcovHAC(model1)),
##     coeftest(model2, vcov=vcovHAC(model2)),
##     report=('vcsp'),
##     type="latex",
##     omit=c("as.character","Constant"),
##     digits=2
##     )
## nobs(model1); nobs(model2)


fromCut<-function(cutoff)
{
	twoallies<-nonUSData[nonUSData[,"idealpointdistance.1"]<=cutoff&nonUSData[,"idealpointdistance.2"]<=cutoff,]
	allyTS<-xts(twoallies[,"GOLDSTEINSCALE"],
              order.by=as.Date(paste(twoallies[,"year"], substr(twoallies$date, 6, 7), "01", sep = "-"))
                    )
	allyMonths<-apply.monthly(allyTS,mean)
        names(allyMonths)[1] <- "conflict"
	aelection<-year(allyMonths)%%4==0
	proximate<-month(allyMonths)%in%c(6,7,8,9,10,11)&aelection
	amidterm<-year(allyMonths)%%4==2
        proximatemidterm <- month(allyMonths)%in%c(6,7,8,9,10,11)&amidterm
        model1<-lm(conflict~proximate+as.character(floor((year(allyMonths)-1)/4))+as.character(month(allyMonths)), data = allyMonths)
        first <- coeftest(model1, vcovHAC(model1))[2,1:2]
        model2<-lm(conflict~proximate+as.character(floor((year(allyMonths)-1)/2))+as.character(month(allyMonths)), data = allyMonths)
        second <- coeftest(model2, vcovHAC(model2))[2,1:2]
	return(c(first,second))
}

## very slow
someCuts<-1.2+0:60/20
res<-mclapply(someCuts,fromCut,mc.cores=2)
res<-do.call(rbind,res)

cutoff<-median(c(nonUSData[,"idealpointdistance.2"],nonUSData[,"idealpointdistance.1"]))

minAlignedEventCoverage <- sapply(c(cutoff, someCuts),function(x) min(tapply(nonUSData[,"idealpointdistance.1"]<=x&nonUSData[,"idealpointdistance.2"]<=x, nonUSData[,"year"], mean)))

maxAlignedEventCoverage <- sapply(c(cutoff, someCuts),function(x) max(tapply(nonUSData[,"idealpointdistance.1"]<=x&nonUSData[,"idealpointdistance.2"]<=x, nonUSData[,"year"], mean)))

meanAlignedEventCoverage <- sapply(c(cutoff, someCuts),function(x) mean(tapply(nonUSData[,"idealpointdistance.1"]<=x&nonUSData[,"idealpointdistance.2"]<=x, nonUSData[,"year"], mean)))


## pdf("cutoffdependence_fouryearblocking.pdf",width=14,height=7)
par(mfrow=c(1,2))

secondCI<-cbind(res[,3]+1.96*res[,4],res[,3]-1.96*res[,4])
plot(minAlignedEventCoverage[-1],res[,3],ylim=c(-0.7,0.3),type="l",main="US-Aligned: 2 year fixed effects, 6 month window", bty="n", xlab="Minimum Proportion of Events in US Bloc in a Year", ylab="Election-Conflict Estimate")
lines(minAlignedEventCoverage[-1],secondCI[,1],col="red",lty=2)
lines(minAlignedEventCoverage[-1],secondCI[,2],col="red",lty=2)
abline(h=0,col="green")
abline(v=minAlignedEventCoverage[1],col="blue")

firstCI<-cbind(res[,1]+1.96*res[,2],res[,1]-1.96*res[,2])
plot(minAlignedEventCoverage[-1],res[,1],ylim=c(-0.7,0.3),type="l",main="US-Aligned: 4 year fixed effects, 6 month window", bty="n", xlab="Minimum Proportion of Events in US Bloc in a Year", ylab="Election-Conflict Estimate")
lines(minAlignedEventCoverage[-1],firstCI[,1],col="red",lty=2)
lines(minAlignedEventCoverage[-1],firstCI[,2],col="red",lty=2)
abline(h=0,col="green")
abline(v=minAlignedEventCoverage[1],col="blue")

## dev.off()


################
cutoff<-median(c(nonUSData[,"idealpointdistance.2"],nonUSData[,"idealpointdistance.1"]))
inUSBloc<-nonUSData[nonUSData[,"ally.1"]&nonUSData[,"ally.2"],]
usBlocTS<-xts(inUSBloc[,"GOLDSTEINSCALE"],
              order.by=as.Date(paste(inUSBloc[,"year"], substr(inUSBloc$date, 6, 7), "01", sep = "-"))
                  )
monthlyMean<-apply.monthly(usBlocTS,mean)
monthlyCount<-apply.monthly(usBlocTS,nrow)
isElection<-year(monthlyCount)%%4==0
inWindow<-month(monthlyCount)%in%c(6,7,8,9,10,11)
monthlyCount<-data.frame(monthlyCount,isElection)

proximateElection<-isElection&inWindow
proximateMidterm<-year(monthlyMean)%%4==2&inWindow
theMonth<-month(monthlyMean)
theYear<-year(monthlyMean)
augmented<-cbind("meanScale"=monthlyMean,"proximate"=proximateElection,"premidterm"=proximateMidterm,"Month"=theMonth,"Year"=theYear)
reelect<-augmented[,"proximate"]&(!augmented[,"Year"]%in%c(1952,1960,1968,1988))
norelect<-augmented[,"proximate"]&(augmented[,"Year"]%in%c(1952,1960,1968,1988))

# Table 19 alliances
model1.a<-lm(meanScale~proximate+as.character(Month),data=augmented,weights=sqrt(monthlyCount$monthlyCount))
coeftest(model1.a,vcov = vcovHAC(model1.a))
#
model2.a<-lm(meanScale~proximate+as.character(Month)+as.character(ceiling(Year/2)),data=augmented,weights=sqrt(monthlyCount$monthlyCount))
coeftest(model2.a,vcov = vcovHAC(model2.a))
#
model3.a<-lm(meanScale~proximate+as.character(Month)+as.character(floor((Year-1)/4)),data=augmented,weights=sqrt(monthlyCount$monthlyCount))
coeftest(model3.a,vcov = vcovHAC(model3.a))

################
cutoff<-median(c(nonUSData[,"idealpointdistance.2"],nonUSData[,"idealpointdistance.1"]))
inUSBloc<-nonUSData[
nonUSData[,"idealpointdistance.1"]<=cutoff&nonUSData[,"idealpointdistance.2"]<=cutoff
               |nonUSData[,"ally.1"]&nonUSData[,"ally.2"]
|nonUSData[,"idealpointdistance.1"]<=cutoff&nonUSData[,"ally.2"]
|nonUSData[,"ally.1"]&nonUSData[,"idealpointdistance.2"]<=cutoff
                 ,]
usBlocTS<-xts(inUSBloc[,"GOLDSTEINSCALE"],
              order.by=as.Date(paste(inUSBloc[,"year"], substr(inUSBloc$date, 6, 7), "01", sep = "-"))
                  )
monthlyMean<-apply.monthly(usBlocTS,mean)
monthlyCount<-apply.monthly(usBlocTS,nrow)
isElection<-year(monthlyCount)%%4==0
inWindow<-month(monthlyCount)%in%c(6,7,8,9,10,11)
monthlyCount<-data.frame(monthlyCount,isElection)

proximateElection<-isElection&inWindow
proximateMidterm<-year(monthlyMean)%%4==2&inWindow
theMonth<-month(monthlyMean)
theYear<-year(monthlyMean)
augmented<-cbind("meanScale"=monthlyMean,"proximate"=proximateElection,"premidterm"=proximateMidterm,"Month"=theMonth,"Year"=theYear)
reelect<-augmented[,"proximate"]&(!augmented[,"Year"]%in%c(1952,1960,1968,1988))
norelect<-augmented[,"proximate"]&(augmented[,"Year"]%in%c(1952,1960,1968,1988))

# Table 19 alliances plus aligned
model1.b<-lm(meanScale~proximate+as.character(Month),data=augmented,weights=sqrt(monthlyCount$monthlyCount))
coeftest(model1.b,vcov = vcovHAC(model1.b))
#
model2.b<-lm(meanScale~proximate+as.character(Month)+as.character(ceiling(Year/2)),data=augmented,weights=sqrt(monthlyCount$monthlyCount))
coeftest(model2.b,vcov = vcovHAC(model2.b))
#
model3.b<-lm(meanScale~proximate+as.character(Month)+as.character(floor((Year-1)/4)),data=augmented,weights=sqrt(monthlyCount$monthlyCount))
coeftest(model3.b,vcov = vcovHAC(model3.b))

## stargazer(
##     coeftest(model1.a, vcov=vcovHAC(model1.a)),
##     coeftest(model2.a, vcov=vcovHAC(model2.a)),
##     coeftest(model3.a, vcov=vcovHAC(model3.a)),
##     coeftest(model1.b, vcov=vcovHAC(model1.b)),
##     coeftest(model2.b, vcov=vcovHAC(model2.b)),
##     coeftest(model3.b, vcov=vcovHAC(model3.b)),
##     report=('vcsp'),
##     type="latex",
##     omit=c("as.character","Constant"),
##     digits=2
##     )
## nobs(model1.a); nobs(model2.a); nobs(model3.a)
## nobs(model1.b); nobs(model2.b); nobs(model3.b)


## COUNTRY JACK-KNIFE

data <- nonUSData

################
cutoff<-median(c(data[,"idealpointdistance.2"],data[,"idealpointdistance.1"]))


inUSBloca<-data[data[,"ally.1"]&data[,"ally.2"],]
inUSBlocb<-data[
data[,"idealpointdistance.1"]<=cutoff&data[,"idealpointdistance.2"]<=cutoff
,]

inUSBlocbAGG1 <- aggregate(
cbind(idealpointdistance.1, 1)~ actor1country,
data = subset(inUSBlocb, idealpointdistance.1 < cutoff),
FUN = sum
)
names(inUSBlocbAGG1)[c(1,3)] <- c("country","n.1")
inUSBlocbAGG2 <- aggregate(
cbind(idealpointdistance.2, 1)~ actor2country,
data = subset(inUSBlocb, idealpointdistance.2 < cutoff),
FUN = sum
)
names(inUSBlocbAGG2)[c(1,3)] <- c("country","n.2")

inUSBlocbAGG <- merge(
inUSBlocbAGG1,
inUSBlocbAGG2
)
inUSBlocbAGG$idealpointdistance <- with(
inUSBlocbAGG,
(idealpointdistance.1 + idealpointdistance.2) / (n.1 + n.2)
)

countryATable <- data.frame(table(c(inUSBloca$actor1country, inUSBloca$actor2country)))
countryATable$Var1 <- as.character(countryATable$Var1)

countryBTable <- data.frame(table(c(inUSBlocb$actor1country, inUSBlocb$actor2country)))
countryBTable$Var1 <- as.character(countryBTable$Var1)

dataDiff <- merge(
countryATable,
countryBTable,
by="Var1",
all=TRUE
)
dataDiff$Freq.x[is.na(dataDiff$Freq.x)] <- 0
dataDiff$diff <- with(dataDiff, Freq.y - Freq.x)
dataDiff <- dataDiff[order(dataDiff$diff,decreasing=T),]
dataDiff$Var1 <- as.character(dataDiff$Var1)
names(dataDiff)[1] <- "country"

dataDiff2 <- merge(
countryATable,
inUSBlocbAGG,
by.x="Var1",
by.y="country",
all=TRUE
)
dataDiff2 <- dataDiff2[order(dataDiff2$idealpointdistance,decreasing=F),]
names(dataDiff2)[1] <- "country"
dataDiff2 <- subset(dataDiff2, country %in% dataDiff$country[1:50])


if (F) {

coefs3 <- NULL

for (i in 1:50) {

inUSBloc<-data[
data[,"idealpointdistance.1"]<=cutoff&data[,"idealpointdistance.2"]<=cutoff
&!(data[,"actor1country"]%in%dataDiff$country[i]&data[,"idealpointdistance.2"]<=cutoff)
&!(data[,"idealpointdistance.1"]<=cutoff&data[,"actor2country"]%in%dataDiff$country[i])
                 ,]

print(dataDiff$country[i])
print(i)
print(nrow(inUSBloc))
usBlocTS<-xts(inUSBloc[,"GOLDSTEINSCALE"],
              order.by=as.Date(paste(inUSBloc[,"year"], substr(inUSBloc$date, 6, 7), "01", sep = "-"))
                  )
monthlyMean<-apply.monthly(usBlocTS,mean)
monthlyCount<-apply.monthly(usBlocTS,nrow)
isElection<-year(monthlyCount)%%4==0
inWindow<-month(monthlyCount)%in%c(6,7,8,9,10,11)
monthlyCount<-data.frame(monthlyCount,isElection)

proximateElection<-isElection&inWindow
proximateMidterm<-year(monthlyMean)%%4==2&inWindow
theMonth<-month(monthlyMean)
theYear<-year(monthlyMean)
augmented<-cbind("meanScale"=monthlyMean,"proximate"=proximateElection,"premidterm"=proximateMidterm,"Month"=theMonth,"Year"=theYear)

#Model 1 of Table 3
coef1.orig<-coef(lm(meanScale~proximate+as.character(Month),data=augmented,weights=sqrt(monthlyCount$monthlyCount)))[2]
se1.orig<-summary(lm(meanScale~proximate+as.character(Month),data=augmented,weights=sqrt(monthlyCount$monthlyCount)))$coefficients[2,2]
#Model 2
coef2.orig<-coef(lm(meanScale~proximate+as.character(Month)+as.character(ceiling(Year/2)),data=augmented,weights=sqrt(monthlyCount$monthlyCount)))[2]
se2.orig<-summary(lm(meanScale~proximate+as.character(Month)+as.character(ceiling(Year/2)),data=augmented,weights=sqrt(monthlyCount$monthlyCount)))$coefficients[2,2]
#Model 3
coef3.orig<-coef(lm(meanScale~proximate+as.character(Month)+as.character(floor((Year-1)/4)),data=augmented,weights=sqrt(monthlyCount$monthlyCount)))[2]
se3.orig<-summary(lm(meanScale~proximate+as.character(Month)+as.character(floor((Year-1)/4)),data=augmented,weights=sqrt(monthlyCount$monthlyCount)))$coefficients[2,2]

print(cbind(coef1.orig, coef2.orig, coef3.orig))

coefs4 <- NULL

n.drops <- nrow(data[
data[,"idealpointdistance.1"]<=cutoff&data[,"idealpointdistance.2"]<=cutoff
                 ,]) - nrow(data[
data[,"idealpointdistance.1"]<=cutoff&data[,"idealpointdistance.2"]<=cutoff
&(data[,"actor1country"]!=dataDiff$country[i]&data[,"actor2country"]!=dataDiff$country[i])
                 ,])

for (j in 1:1000) {

inUSBloc<-data[
data[,"idealpointdistance.1"]<=cutoff&data[,"idealpointdistance.2"]<=cutoff
                 ,]
inUSBloc<- inUSBloc[sample(nrow(inUSBloc),size=nrow(inUSBloc)-n.drops),]

usBlocTS<-xts(inUSBloc[,"GOLDSTEINSCALE"],
              order.by=as.Date(paste(inUSBloc[,"year"], substr(inUSBloc$date, 6, 7), "01", sep = "-"))
                  )

monthlyMean<-apply.monthly(usBlocTS,mean)
monthlyCount<-apply.monthly(usBlocTS,nrow)
isElection<-year(monthlyCount)%%4==0
inWindow<-month(monthlyCount)%in%c(6,7,8,9,10,11)
monthlyCount<-data.frame(monthlyCount,isElection)

proximateElection<-isElection&inWindow
proximateMidterm<-year(monthlyMean)%%4==2&inWindow
theMonth<-month(monthlyMean)
theYear<-year(monthlyMean)
augmented<-cbind("meanScale"=monthlyMean,"proximate"=proximateElection,"premidterm"=proximateMidterm,"Month"=theMonth,"Year"=theYear)

#Model 1 of Table 3
coef1<-coef(lm(meanScale~proximate+as.character(Month),data=augmented,weights=sqrt(monthlyCount$monthlyCount)))[2]
#Model 2
coef2<-coef(lm(meanScale~proximate+as.character(Month)+as.character(ceiling(Year/2)),data=augmented,weights=sqrt(monthlyCount$monthlyCount)))[2]
#Model 3
coef3<-coef(lm(meanScale~proximate+as.character(Month)+as.character(floor((Year-1)/4)),data=augmented,weights=sqrt(monthlyCount$monthlyCount)))[2]

coefs4<- rbind(
    coefs4,
    data.frame(
        dataDiff$country[j],
        j,
        nrow(inUSBloc),
        coef1, coef2, coef3
        )
    )

if (j%%100==0) {
    print(j)
}

           }

coefs3 <- rbind(
    coefs3,
    data.frame(
        dataDiff$country[i],
        i,
        nrow(inUSBloc)-n.drops,
        coef1.orig, coef2.orig, coef3.orig, se1.orig, se2.orig, se3.orig,
        lower1=quantile(coefs4$coef1, probs=0.025),
        upper1=quantile(coefs4$coef1, probs=0.975),
        lower2=quantile(coefs4$coef2, probs=0.025),
        upper2=quantile(coefs4$coef2, probs=0.975),
        lower3=quantile(coefs4$coef3, probs=0.025),
        upper3=quantile(coefs4$coef3, probs=0.975)
        )
    )

save(coefs3,file="country_removal_significant.RData")

pdf("alliance_expand_by_country.pdf",width=14)
par(mfrow=c(1,2))

plot(coefs3$coef2,ylim=c(-0.7,0.2),main="US-Allied (x=0) + US-Aligned, by N events added\n2 year fixed effects, 6 month window", bty="n", xlab="From Alliances Only (x=0), Add Nth Country \n(most added events to least added events)", ylab="Election-Conflict Estimate", pch=16)
segments(x0=1:nrow(coefs3),y0=coefs3$lower2,x1=1:nrow(coefs3),y1=coefs3$upper2)
abline(h=0,col="green")
text(1:nrow(coefs3),coefs3$coef2,labels=dataDiff$country[1:nrow(coefs3)],cex=0.5,pos=3)

plot(coefs3$coef3,ylim=c(-0.7,0.2),main="US-Allied (x=0) + US-Aligned, by N events added\n4 year fixed effects, 6 month window", bty="n", xlab="From Alliances Only (x=0), Add Nth Country \n(most added events to least added events)", ylab="Election-Conflict Estimate", pch=16)
segments(x0=1:nrow(coefs3),y0=coefs3$lower3,x1=1:nrow(coefs3),y1=coefs3$upper3)
abline(h=0,col="green")
text(1:nrow(coefs3),coefs3$coef3,labels=dataDiff$country[1:nrow(coefs3)],cex=0.5,pos=3)

dev.off()

           }

}

## YEAR JACK-KNIFE
################
cutoff<-median(c(data[,"idealpointdistance.2"],data[,"idealpointdistance.1"]))


post.election <- unique(floor((1946:1993-1)/4)*4)+1

THE.YEARS <- lapply(post.election, function(x) x:(x+3))

if (F) {

coefs3 <- NULL

for (i in THE.YEARS) {

inUSBloc<-data[
data[,"idealpointdistance.1"]<=cutoff&data[,"idealpointdistance.2"]<=cutoff
&!(data[,"year"]%in%i)
                 ,]

print(i)
print(nrow(inUSBloc))
usBlocTS<-xts(inUSBloc[,"GOLDSTEINSCALE"],
              order.by=as.Date(paste(inUSBloc[,"year"], substr(inUSBloc$date, 6, 7), "01", sep = "-"))
                  )
monthlyMean<-apply.monthly(usBlocTS,mean)
monthlyCount<-apply.monthly(usBlocTS,nrow)
isElection<-year(monthlyCount)%%4==0
inWindow<-month(monthlyCount)%in%c(6,7,8,9,10,11)
monthlyCount<-data.frame(monthlyCount,isElection)

proximateElection<-isElection&inWindow
proximateMidterm<-year(monthlyMean)%%4==2&inWindow
theMonth<-month(monthlyMean)
theYear<-year(monthlyMean)
augmented<-cbind("meanScale"=monthlyMean,"proximate"=proximateElection,"premidterm"=proximateMidterm,"Month"=theMonth,"Year"=theYear)

#Model 1 of Table 3
coef1.orig<-coef(lm(meanScale~proximate+as.character(Month),data=augmented,weights=sqrt(monthlyCount$monthlyCount)))[2]
se1.orig<-summary(lm(meanScale~proximate+as.character(Month),data=augmented,weights=sqrt(monthlyCount$monthlyCount)))$coefficients[2,2]
#Model 2
coef2.orig<-coef(lm(meanScale~proximate+as.character(Month)+as.character(ceiling(Year/2)),data=augmented,weights=sqrt(monthlyCount$monthlyCount)))[2]
se2.orig<-summary(lm(meanScale~proximate+as.character(Month)+as.character(ceiling(Year/2)),data=augmented,weights=sqrt(monthlyCount$monthlyCount)))$coefficients[2,2]
#Model 3
coef3.orig<-coef(lm(meanScale~proximate+as.character(Month)+as.character(floor((Year-1)/4)),data=augmented,weights=sqrt(monthlyCount$monthlyCount)))[2]
se3.orig<-summary(lm(meanScale~proximate+as.character(Month)+as.character(floor((Year-1)/4)),data=augmented,weights=sqrt(monthlyCount$monthlyCount)))$coefficients[2,2]

print(cbind(coef1.orig, coef2.orig, coef3.orig))

coefs4 <- NULL

n.drops <- nrow(data[
data[,"idealpointdistance.1"]<=cutoff&data[,"idealpointdistance.2"]<=cutoff
                 ,]) - nrow(data[
data[,"idealpointdistance.1"]<=cutoff&data[,"idealpointdistance.2"]<=cutoff
&!(data[,"year"]%in%i)
                 ,])

for (j in 1:1000) {

inUSBloc<-data[
data[,"idealpointdistance.1"]<=cutoff&data[,"idealpointdistance.2"]<=cutoff
                 ,]
inUSBloc<- inUSBloc[sample(nrow(inUSBloc),size=nrow(inUSBloc)-n.drops),]

## print(dataDiff$country[j])
## print(i)
## print(nrow(inUSBloc))

usBlocTS<-xts(inUSBloc[,"GOLDSTEINSCALE"],
              order.by=as.Date(paste(inUSBloc[,"year"], substr(inUSBloc$date, 6, 7), "01", sep = "-"))
                  )

monthlyMean<-apply.monthly(usBlocTS,mean)
monthlyCount<-apply.monthly(usBlocTS,nrow)
isElection<-year(monthlyCount)%%4==0
inWindow<-month(monthlyCount)%in%c(6,7,8,9,10,11)
monthlyCount<-data.frame(monthlyCount,isElection)

proximateElection<-isElection&inWindow
proximateMidterm<-year(monthlyMean)%%4==2&inWindow
theMonth<-month(monthlyMean)
theYear<-year(monthlyMean)
augmented<-cbind("meanScale"=monthlyMean,"proximate"=proximateElection,"premidterm"=proximateMidterm,"Month"=theMonth,"Year"=theYear)

#Model 1 of Table 3
coef1<-coef(lm(meanScale~proximate+as.character(Month),data=augmented,weights=sqrt(monthlyCount$monthlyCount)))[2]
#Model 2
coef2<-coef(lm(meanScale~proximate+as.character(Month)+as.character(ceiling(Year/2)),data=augmented,weights=sqrt(monthlyCount$monthlyCount)))[2]
#Model 3
coef3<-coef(lm(meanScale~proximate+as.character(Month)+as.character(floor((Year-1)/4)),data=augmented,weights=sqrt(monthlyCount$monthlyCount)))[2]

## print(cbind(coef1, coef2, coef3))

coefs4<- rbind(
    coefs4,
    data.frame(
        i[4],
        j,
        nrow(inUSBloc),
        coef1, coef2, coef3
        )
    )

if (j%%100==0) {
    print(j)
}

           }

coefs3 <- rbind(
    coefs3,
    data.frame(
        i[4],
        nrow(inUSBloc)-n.drops,
        coef1.orig, coef2.orig, coef3.orig, se1.orig, se2.orig, se3.orig,
        lower1=quantile(coefs4$coef1, probs=0.025),
        upper1=quantile(coefs4$coef1, probs=0.975),
        lower2=quantile(coefs4$coef2, probs=0.025),
        upper2=quantile(coefs4$coef2, probs=0.975),
        lower3=quantile(coefs4$coef3, probs=0.025),
        upper3=quantile(coefs4$coef3, probs=0.975)
        )
    )

save(coefs3,file="year_removal_significant.RData")

}}


## Presidential approval

## "presidential_approval.txt" is Gallup survey data
## DATA <- read.table("presidential_approval.txt", head=FALSE)
## for_fill is a placeholder for the same president as the previous period
## DATA[DATA$V1=="for_fill","V1"] <- NA
## names(DATA) <- c(
##     "president","start_date","end_date",
##     "approving","disapproving","unsure_no_data"
##     )

## DATA$president <- na.locf(DATA$president)

## DATA$start_date <- as.Date(DATA$start_date, format="%m/%d/%Y")
## DATA$start_year <- year(DATA$start_date)
## DATA$start_month <- month(DATA$start_date)
## DATA$start_date2 <- with(DATA, paste(start_year, start_month, "01", sep = "-"))

## AGG <- aggregate(approving ~ start_date2, data = DATA, FUN = mean)

## TS <- with(AGG, zoo(approving, as.Date(start_date2, format="%Y-%m-%d")))
## TS <- ## na.approx(
##     merge(
##         TS,
##         zoo(,seq(as.Date("1941-07-01"), as.Date("1993-12-31"), by = "month"))
##         )
##     ## )


## cutoff<-median(c(nonUSData[,"idealpointdistance.2"],nonUSData[,"idealpointdistance.1"]))
## inUSBloc<-nonUSData[nonUSData[,"idealpointdistance.1"]<=cutoff&nonUSData[,"idealpointdistance.2"]<=cutoff,]
## #
## usBlocTS<-xts(inUSBloc[,"GOLDSTEINSCALE"],
##               order.by=as.Date(paste(inUSBloc[,"year"], substr(inUSBloc$date, 6, 7), "01", sep = "-"))
##                   )
## monthlyMean<-apply.monthly(usBlocTS,mean)
## monthlyCount<-apply.monthly(usBlocTS,nrow)
## isElection<-year(monthlyCount)%%4==0
## inWindow<-month(monthlyCount)%in%c(6,7,8,9,10,11)
## proximateElection<-isElection&inWindow
## proximateMidterm<-year(monthlyMean)%%4==2&inWindow
## theMonth<-month(monthlyMean)
## theYear<-year(monthlyMean)
## augmented<-cbind("meanScale"=monthlyMean,"proximate"=proximateElection,"premidterm"=proximateMidterm,"Month"=theMonth,"Year"=theYear)

## approvalAndConflict <- merge(
##     TS[index(TS) >= "1946-01-01" & index(TS) <= "1993-12-31"],
##     monthlyMean[index(monthlyMean) >= "1946-01-01" & index(monthlyMean) <= "1993-12-31",]
##     )
## approvalAndConflict <- merge(
##     approvalAndConflict,
##     monthlyCount[index(monthlyCount) >= "1946-01-01" & index(monthlyCount) <= "1993-12-31",]
##     )
## names(approvalAndConflict) <- c("approval", "conflict", "conflictcount")

## isElection<-year(approvalAndConflict)%%4==0
## inWindow<-month(approvalAndConflict)%in%c(6,7,8,9,10,11)

## proximateElection<-isElection&inWindow
## proximateMidterm<-year(approvalAndConflict)%%4==2&inWindow
## theMonth<-month(approvalAndConflict)
## theYear<-year(approvalAndConflict)
## lagappr<-c(NA,approvalAndConflict[1:(nrow(approvalAndConflict)-1),1])
## augmented<-cbind("meanScale"=approvalAndConflict$conflict,"confCount"=approvalAndConflict$conflictcount,"approval"=approvalAndConflict$approval,"proximate"=proximateElection,"premidterm"=proximateMidterm,"Month"=theMonth,"Year"=theYear)
## augmented$reelect<-augmented[,"proximate"]&(!augmented[,"Year"]%in%c(1952,1960,1968,1988))
## augmented$norelect<-augmented[,"proximate"]&(augmented[,"Year"]%in%c(1952,1960,1968,1988))

## model1<-lm(meanScale~proximate*scale(approval)+as.character(Month),data=augmented)
## coeftest(model1,vcov = vcovHAC(model1))

## stargazer(
##     coeftest(model1, vcov=vcovHAC(model1)),
##     report=('vcsp'),
##     type="latex",
##     omit=c("as.character","Constant"),
##     digits=2
##     )
## nobs(model1)
